library(rgdal)
library(leaflet)
library(sp)
library(dplyr)
#Unzip files
unzip("CO-wells.zip", exdir = ".")
list.files("CO-wells")
## [1] "CO-HUC.cpg"   "CO-HUC.dbf"   "CO-HUC.prj"   "CO-HUC.qpj"  
## [5] "CO-HUC.shp"   "CO-HUC.shx"   "CO-wells.csv"
#Read in data and convert to SpatialPointsDataFrame
CO.Wells <- read.csv(file ="CO-wells/CO-wells.csv", header = T, na.strings = "NA")
wells <- SpatialPointsDataFrame(coords = CO.Wells[,2:1], data = CO.Wells, proj4string = CRS("+proj=longlat +datum=WGS84"))

#Read in polygons and change change proj4string
huc1 <- readOGR("CO-wells", "CO-HUC")
## OGR data source with driver: ESRI Shapefile 
## Source: "CO-wells", layer: "CO-HUC"
## with 93 features
## It has 2 fields
huc <- spTransform(huc1, CRS("+proj=longlat +datum=WGS84"))

#Records HUC ID for each point in wells
wells@data$HUC_ID <- over(wells, huc)$HUC_ID

#Groups by HUC ID and records total wells, active, and inactive
huc.wells <- wells@data %>% group_by(HUC_ID) %>% summarize(total = length(Facil_Stat), active = sum(Facil_Stat == "AC" | Facil_Stat == "PR"), inactive = total  - active)
## Warning: package 'bindrcpp' was built under R version 3.4.2
huc.wells
## # A tibble: 91 x 4
##    HUC_ID total active inactive
##     <int> <int>  <int>    <int>
##  1    873     1      0        1
##  2    880     5      0        5
##  3    899   773    179      594
##  4    903   313    139      174
##  5    932   203     31      172
##  6    936     2      0        2
##  7    947  3095   1404     1691
##  8    951  4308   1945     2363
##  9    963  1231    571      660
## 10    964     4      0        4
## # ... with 81 more rows
#Joins huc.wells and the data of huc by HUC_ID
huc@data <- left_join(huc@data, huc.wells, by = "HUC_ID")
dim(huc@data)
## [1] 93  5
head(huc@data)
##   HUC_ID                            HUC_NAME total active inactive
## 1    846                  Upper North Platte    NA     NA       NA
## 2    873 Upper Green-Flaming Gorge Reservoir     1      0        1
## 3    880                       Upper Laramie     5      0        5
## 4    899                        Little Snake   773    179      594
## 5    903                           Vermilion   313    139      174
## 6    932                     Lower Lodgepole   203     31      172
#Replace NA values with 0
huc@data[is.na(huc@data)] <- 0
huc@data
##    HUC_ID                            HUC_NAME total active inactive
## 1     846                  Upper North Platte     0      0        0
## 2     873 Upper Green-Flaming Gorge Reservoir     1      0        1
## 3     880                       Upper Laramie     5      0        5
## 4     899                        Little Snake   773    179      594
## 5     903                           Vermilion   313    139      174
## 6     932                     Lower Lodgepole   203     31      172
## 7     936                     Upper Lodgepole     2      0        2
## 8     947                                Crow  3095   1404     1691
## 9     951                     Cache La Poudre  4308   1945     2363
## 10    963                       Lone Tree-Owl  1231    571      660
## 11    964                 Lower Green-Diamond     4      0        4
## 12    972                  Lower South Platte    35      0       35
## 13    977                         Sidney Draw   281      5      276
## 14    983             North Platte Headwaters   756     98      658
## 15    985                         Upper Yampa   938    122      816
## 16    989                         Lower Yampa   591     45      546
## 17    996        Middle South Platte-Sterling  6536    272     6264
## 18   1007                              Pawnee  2970    646     2324
## 19   1011                      Stinking Water   171     69      102
## 20   1022                         Lower White  4601   1354     3247
## 21   1032                         Upper White   522    134      388
## 22   1033                        Big Thompson  2228   1142     1086
## 23   1034                    Upper Republican     1      0        1
## 24   1036    Middle South Platte-Cherry Creek 23080  12248    10832
## 25   1038                           Frenchman   391    176      215
## 26   1051                 Colorado headwaters    77      0       77
## 27   1067                     Piceance-Yellow  2239    795     1444
## 28   1068               North Fork Republican  6375   2682     3693
## 29   1086                           St. Vrain  4944   2868     2076
## 30   1097                               Kiowa  1188    197      991
## 31   1100                              Beaver  3179    110     3069
## 32   1105                               Bijou  1809    121     1688
## 33   1114                                Blue     0      0        0
## 34   1124         Colorado headwaters-Plateau 12615   8328     4287
## 35   1128                      Parachute-Roan  5575   3099     2476
## 36   1131                               Eagle     2      0        2
## 37   1138                               Clear    45      3       42
## 38   1141                            Arikaree  2335   1086     1249
## 39   1149               South Fork Republican   800    313      487
## 40   1150                    Westwater Canyon    10      0       10
## 41   1152                  Upper South Platte    24      0       24
## 42   2207                        Roaring Fork    33     10       23
## 43   1161                       Little Beaver    40      0       40
## 44   1173                   South Fork Beaver    71      1       70
## 45   1176             South Platte Headwaters    28      0       28
## 46   1182                 Arkansas Headwaters    17      0       17
## 47   1187                 North Fork Gunnison   187     30      157
## 48   1196                      Lower Gunnison   222      8      214
## 49   1207                           Big Sandy  1560    210     1350
## 50   1215                       Lower Dolores     2      0        2
## 51   2208                         East-Taylor     1      0        1
## 52   1225                            Fountain    49      0       49
## 53   1226                                Rush   355     27      328
## 54   1234               North Fork Smoky Hill   111      2      109
## 55   1236                               Chico    56      0       56
## 56   1240                      Upper Gunnison    12      0       12
## 57   1241                        Uncompahange    48      0       48
## 58   1245                               Horse    82      0       82
## 59   1250               Smoky Hill Headwaters   493    100      393
## 60   1258                      Upper Arkansas   424     59      365
## 61   1264                       Upper Dolores   205     10      195
## 62   1265                          San Miguel   246     86      160
## 63   1266                             Tomichi     1      0        1
## 64   1277        Upper Arkansas-Lake Meredith    51      0       51
## 65   1278                              Ladder   440    100      340
## 66   1282          Upper Arkansas-John Martin  1053     86      967
## 67   1291                          Whitewoman   239     19      220
## 68   1295                            San Luis    22      0       22
## 69   1302                            Saguache     3      0        3
## 70   1308                           Montezuma   211     31      180
## 71   1319         Lower San Juan-Four Corners   136      4      132
## 72   1328       Middle Arkansas-Lake Mckinney     8      0        8
## 73   1331                            Huerfano   371     37      334
## 74   1332                              Animas  1765   1294      471
## 75   1336               Rio Grande headwaters     4      0        4
## 76   1341                            Apishapa   426    241      185
## 77   1350                   Alamosa-Trinchera    16      0       16
## 78   1352                          Purgatoire  3531   2516     1015
## 79   1353                              Mcelmo   366     95      271
## 80   1355                      Upper San Juan  2452   1547      905
## 81   1361                           Two Butte    64      1       63
## 82   1365                              Piedra    90     38       52
## 83   1371                              Mancos   463     10      453
## 84   1380                                Bear   137      3      134
## 85   1381                     Middle San Juan   761    338      423
## 86   1400                             Conejos     4      0        4
## 87   1402                 North Fork Cimarron   335     46      289
## 88   1415                         Sand Arroyo   167     36      131
## 89   1423                           Rio Chama     3      0        3
## 90   1426                      Upper Cimarron   307     75      232
## 91   1432                    Upper Rio Grande     2      0        2
## 92   1435                 Cimarron headwaters    13      0       13
## 93   1444                 Canadian headwaters    19      0       19
pal <- colorBin("YlOrRd", domain = huc$total, bins = 6, pretty = TRUE)
#Create custom labels
labels <- sprintf(
  "<strong>%s</strong><br/>%s<br/>%s<br/>%s",
  paste("HUC:",huc@data$HUC_NAME), paste("Total:",huc@data$total), paste("Active:",huc@data$active), paste("Inactive:",huc@data$inactive)
) %>% lapply(htmltools::HTML)

leaflet(width = "100%") %>% 
# Overlay groups
addPolygons(data = huc, fillColor = ~pal(total), group = "Oil & gas wells", fillOpacity = .4, color = "white", opacity = 1, weight = 2, dashArray = "3", highlight = highlightOptions(
  weight = 5,
  color = "#666",
  dashArray = "",
  fillOpacity = 0.4,
  bringToFront = TRUE),
label = labels,
labelOptions = labelOptions(
  style = list("font-weight" = "normal", padding = "3px 8px"),
  textsize = "15px",
  direction = "auto",
  noHide = TRUE))%>% 
    
# Base groups
addProviderTiles(providers$OpenStreetMap.Mapnik, group = "OSM (default)") %>% 
    addProviderTiles(providers$Esri.WorldImagery, group = "Imagery") %>% addProviderTiles(providers$Esri.WorldTopoMap, 
    group = "Topo") %>% 
# Layers control
addLayersControl(overlayGroups = c("Oil & gas wells"), baseGroups = c("OSM (default)", 
    "Imagery", "Topo"), options = layersControlOptions(collapsed = FALSE)) %>%
#Add Legend
addLegend("bottomright", pal = pal, values = huc$total, title = "Total number of gas and oil<br>wells by USGS HUC")